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The  Goal 


To  consider  the  Magnus  series  expansion  for 
numerically  solving  the  non-linear  ODE 


-7-  =  f(x,£),  x(0)  =  xq  G 


The  Magnus  expansion  and  other  geometric 
numerical  integrators  have  been  successfully 
used  to  solve  the  following  eqs.  separately 


Different  versions  of  the  ME  have  proved 
highly  efficient  in  a  number  of  problems: 

-  Quantum  Mechanics 

-  Highly  oscillatory  problems 

-  Eigenvalue  problems 

-  etc. 


Is  it  possible  to  use  Magnus  for  the  non- 
autonomous  non-linear  problem?  (Zanna) 

Is  it  possible  to  combine  both  of  them  for 
numerically  solving  the  problem? 
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Standard  Procedure:  to  consider  t  =  xt  on 
the  vector-field  and  to  solve  the  autonomous 
equation 


dy 

dt 


=  F(y) 


d 

dt 


using  an  standard  algorithm 


f(x,xt) 
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Problems: 

(i)  In  general,  many  evaluations  of  f(x,t)  at 

different  times,  The  algorithm  can 

be  expensive. 

(ii)  If  the  t-dependent  functions  in  f(x,t)  and 
the  x-dependent  functions  in  f(x,t)  evolve  at 
different  time-scales  it  seems  convenient 
to  treat  them  differently. 

(iii)  Some  times,  the  structure  of  the  vector- 
field  f(x,t)  is  simpler  than  the  structure  of 
F(y),  allowing  to  use  more  efficient  algorithms 

Solution?  To  use  explicit  or  implicit  RK  meth¬ 
ods  with  Gaussian  or  Lobatto  quadrature  points. 
This  solves  (i)  but  not  (ii)  and  (iii) 
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Expected  Advantages  on  using  Magnus 


(i)  Efficient  treatment  of  the  explicit  time- 
dependent  functions  on  the  vector-field  f(x,t) 
in  the  following  cases: 

*  cost(f(x,t!)+f(x,t2))  >  cost(f(x1,t)+f(x2,t)) 

i-e.  x  =  fo(t)  +  +  f2(t)x2  +  f3  (t)x3 

Given  with  ^  —  to  +  Cih,  and 

the  quadrature  points  at  order  n  =>  with 
f(x,ti), . . . ,  f(x,ts)  allows  methods  of  order  n 

(ii)  The  f-dependent  functions  in  f(x,t)  and 
the  x-dependent  functions  in  f(x,t)  are  treat¬ 
ed  separately. 

(iii)  It  is  of  particular  interest  for  those  prob¬ 
lems  where  efficient  integrators  for  the  au¬ 
tonomous  vector-field  f(x,r)  are  known  (with 
r  a  constant). 

(iv)  The  geometric  properties  are  preserved. 
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Magnus  expansion  for  linear  systems 


dx 

dt 


—  A(t)x 


x(/i)  =  e  xq, 


where  Q  =  ^/,  with 


» h 


Qi  =  /  A(t)dt 


o 


\  [  dti  [  dt2[Atl,At2] 

2  Jo  Jo 

—  J  dt\  J  dt2  J  Jt3|  [A^,  At2,  At3]  +  [At3,  At2,  j 

ph  rti  pt2  rt3  , 

— -  /  dti  dt2  /  Jt3  /  [At1?  At4,  At3,  At2] 

12  Jo  Jo  Jo  Jo  ^ 

+  [T2,  T35  At4,  A^]  +  [Atl,  At2 ,  A*3,  TJ  —  [At4,  At3,  T2,  At J  | 


^2 

^3 

^4 


where  Ati  =  A(ti),  [Ai,a2,a3]  =  [Ai,[a2,a3]]  ,  etc. 


This  expansion  is  of  great  interest.  However, 
for  building  efficient  integrators  it  is  conve¬ 
nient  to  write  the  truncated  Magnus  series  in 
a  more  appropriate  way 
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To  take  profit  of  the  time-symmetry,  let  us 
consider  the  Taylor  expansion  of  A(t)  around 
f  1/2  —  to  +  h/2, 


oo 


Mi)  =  ai  / 

i=0 


1  dM(f) 
*!  df 


* — *1/2 


then,  taking  6j  =  i  =  1, 2, 3 


Q 

Q 


6i-— [6i,62]  +  oo5) 

bi  +  T3  -  T[6i,  ft2]  +  T[62)  53] 

1  [6i,  [61,63]]  - -/—[62,  [6i,62]] 


360 

1 

720 


240 


[61,[&i,[61,62]]]+0(/i7) 


{61,62,63}  generators  of  a  graded  free  Lie  al¬ 


gebra  with  grades  1,  2,  3  (Munthe-Kaas  &  Owren). 
[6*,6j,6fc]  =  [6i,[6j,6fe]],  of  order  0{hl+3+k)- 


It  is  easier  to  work  with. 

The  time-symmetry  is  clear 
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In  practice,  it  seems  convenient  to  consider 


A i  =  A(t\),  ...  ,  As  =  A{tg) 
at  some  quadrature  points  with 


Then  3  constants  (3K-  J  such  that 


i  —  0, 1,2  (iserles  &l  N0rsett;  B,  Casas  &l  Ros  ) 
Up  to  order  6:  e  =e  +  0(/in+1) 

nig  =  i 

flW  =  A<°>  +  "I 

Q  =  -  6A&>] 


q[6]  =  #)+Q+[A(0),[A(°),3(2)-  1.  Q]]  +  |[A(1),Q] 

2  60  5 
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This  is  valid  for  any  n-th  order  quadrature. 
We  can  also  use  different  quadrature  formu¬ 
las  for  each  component,  (t)  ■ 

Examples:  Gaussian  Quadratures 

Fourth  order  Consider  At  —  A(cih),  i  —  1,2 
with  c\  2  —  \  T  Then 

=  —  (Ai  +  A2 )  ~  j  A(t)dt 

aW  =  ^r(A2-Ai)*lL 


Sixth  order  Consider  Ai  —  A(q/r),  i  —  1,2,3 
with  C13  =  j  C2  =  1/2.  Then 


A(°) 


h  fh 

(5Ai  T  8A2  T  5A3)  ~  /  A(t)dt 

Jo 


18 
VT5/1 


36 


(A3  —  Ai) 


h 

24 


(Ai  +  A3) 


h 

1  [h 

hj0 


U  H) 


A(t)dt 


t  — 


h 


A  {t)dt 
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This  is  equivalent  to: 


1.  Consider  the  equation:  x7  =  f(x,  t)=  A(t)x 

2.  Frozen  the  x  coordinates  of  the  vector- 
field,  f(x,£) 

3.  Take  the  Magnus  time-average 
fQ(x)=  Q(/i)x  (autonomous  vector-field) 

4.  Consider  the  equation:  y'  =  fQ(y)=  Q(/i)y 

5.  Evaluate  the  1-flow:  y(l)  =  eQ  x(/i) 


Question:  Is  it  possible  to  follow  these  steps 
in  a  simple  way  for  the  non-linear  problem? 
Answer:  YES,  but  only  for  the  steps  1  to  4. 
In  general,  due  to  the  complexity  of  fQ,  step 
5  requires  the  evaluation  of  a  very  complicat¬ 
ed  map 

Solution:  To  approximate  the  complicate  map 
by  a  composition  of  simpler  maps 
Composition  Magnus  Integrators 
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Non-linear  autonomous  system 


x  =  f(x)  (1) 

with  solution  x(t)  =  Ot(xo).  Let  us  consider 
the  Lie  operador 


Li  =  E  fi 

i—  1 

Eg.  (i)  can  be  written  as 


x  =  Lf(x)x 


d 

dt 


<£f(xo)  =  £f<*>f(xo) 


and  in  terms  of  the  evolution  operator 


-  0fLf(xo) 


We  have  rewritten  eq.  (1)  as  a  (infinite  di¬ 
mensional)  linear  eq.  with  formal  solution: 
Oj,  =  exp(tLf) 

Sol.  dif.  eq.  =  Lie  Transformation 


x(t)  =  exp(iLf(y))(y) 


y=xo 
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Non-autonomous  System 


dx. 

dt 


=  f(x,£) 


with  solution  x(t)  =  0^(xq).  The  evolution 
operator  satisfy  the  eq. 


It  is  possible  to  make  use  of  the  Magnus  ex¬ 
pansion  and  to  write  the  solution  as 

0f  =  exp(%>(x0.t))’  with  fQ  =  E  fP 

i 

where 


fp(xo  ,t) 

fp(x0,t) 
being  h  = 


=  /  f  (X0;  s)d-S 

1  ft  fS  1 

=  — —  J  ds\  J  ds2(f(x0,s1),f(x0,S2)) 
:  (f,  g)  the  Lie  bracket 


n  / 

hi  —  (f,  g  )i  —  Lfgi-Lgfi  =  [fj 

j= 1  ' 


dgj 

dxj 


dh 

dxj 
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The  formal  solution  at  t  =  h  is  given  by 


x(/i)  =  exp(Lftlf  ,0(y) 

y=X0 

the  1-flow  solution  of  the  autonomous  DE 


X  =  fQ(x,  h), 


x(0)  =  x0 


f(0  = 

3  = 1 

i  —  0, 1,2 


hl 


'  t--\  f (x,t)dt  +  0(hn+1). 


'0 


Up  to  order  6:  exp  (l{  j  =  exp  (Lfn)+C>(fon"*"1) 


fn[2]  _  f(o) 


f^[4]  —  f(0)  _  (f(l)  f(°)) 


Q  =  -(f(D,|f(o)  _  6f<2>) 


fn[6]  =  f(0)+Q+(f(0);(f(0);if(2)_  J_ Q))  +  |(fd),Q) 

2  60  5 
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The  vector-field  containing  Lie  brackets  are, 
in  general,  very  complicate  and  their  maps 
are  computationally  expensive. 

Solution:  To  approximate  this  map  by  a  com¬ 
position  of  simpler  maps 

This  is  closely  related  to  the  commutator- 
free  methods  (Celledoni  &  Owren) 

Which  simpler  maps  can  be  used? 

A  map  is  the  1-flow  solution  of  a  differential 
equation. 

Which  differential  equations  can  be  exactly 
or  efficiently  computed? 

This  dependes  on  the  particular  problem. 
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Given  the  quadrature  points,  t\, . . . ,  ts,  we  con¬ 
sider  the  following  cases: 

(i)  Suppose  the  autonomous  equation  can  be 
efficiently  solved  up  to  t  —  h 

x  =  aif(x,ti)  -| - 1-  asf(x,ts) 

i.e.  suppose  it  is  easy  to  approximate  the  map 

x(fi)  =  exp  (/iaiZ'f(x,£1)  +  ‘  •  --|-fiasLf(X)ts^(xo) 

(ii)  Given  f  =  +  fg,  suppose  we  can  solve 

x  =  - fasf4(x,  ts) 

x  =  /3ifB(x,ti)4 - f/3sfB(x,fs) 

i.e.  suppose  we  can  approximate  the  maps 

x(/i)  =  exp  (haiLf^j^  -\ - f-  hasLiA^x  tgyj  (x0) 

x(/r)  =  exp  4 - h  /r/3sLfB(x  ts)j(xo) 
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Problem  to  solve 


exp(Q) 


exp(Q) 


exp  (&i  -  ^[61,  M)  +  0(h5) 

m 

exp  (0:^161  +  0^,262)  +  0(h 5) 


i=  1 


exp  ^61  +  ^_63  -  ^[61, 62]  +  ^^[62, 63] 


360 

1 


12 

[61,  [61,63]] 


12 


240 


240 
[62,  [61,62]] 


720 


[61,  [61,  [61,  62]]]  )  +  0(h7) 


m 


[  exp  (ai,i6i  +  0^,262  +  0^,363)  +  0(h7) 


i— 1 


The  coefficients  ot^j  have  to  solve  a  system 
of  non-linear  equations. 


Time-symmetry.  It  can  be  preserved  with 

^771+1—1,1  ^2,1 

^m+l—  i,2  °^,2 

^m+l— i,3  ^i,3 
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Order  Conditions 


Given  the  time-symmetric  composition 

e«l,l&l+al, 2^2+^1,363  e«2,l6l+a2, 262+^2, 363  .  .  . 

.  .  .  e«2,l6l— «2, 262+«2, 363  e«l,l^l— «1, 2&2+<*l, 363 

and  a  general  (time-symmetric)  element 

=  p[k)b!  +  /4fe) b3  +  (3^k)[  1, 2]  +  /3<fe)[2,  3] 

+  4/c)[l,  1,3]  +  /$°[2, 1,2]+  4k)[l,  1,1,2] 

the  order  conditions  can  be  easily  obtained 
from  the  recursive  relation 

e®&!+2/&2+2&3  exbi~yb2+zb3  _  eC(/5^+1) (x,y,z)) 


Next,  we  must  solve  the  equations 

(Mm) ,  4m) 4m) ,  /?im) ,  /4m) ,  4m) >  4m) ) 


111 
12’  240’ 360’ 


1  1  ^ 

240’  720/ 


Finally,  we  write  the  solution  in  terms  of 

A(0)?A(!)?A(2) 
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Example:  Fourth-order  composition  methods 
(B  &  Moan) 


Linear  problem 

en  ~  exp  (bx  +  j^lb2,  &l]) 

~  exp  (a(°)  +  [a(i),a(0)]) 

~  exp  _|_  2A(1A  exp  (L4O)  -  2 

~  exp  (A(1))  exp  (a(°A  exp  (—  aW) 

Non-linear  problem 

exp(Lfn)  ~  exp  (£f(o)_(f(i)f(o))) 

-  exP  (Lif(0)_2f(l))  exp  (Lif(0)+2f(i)) 
^  exp  (— Lf(i))  exp  (Lf(0))  exp  (Lf(i)) 
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X  =  f(x,  t) 


x2 ( 1 )  —  exP  f-^lf(0)_2f(l)J  exp  f-^lf(0)_|_2f(1))  (Xp) 


XI  =  ^f(0)(x!)  -  2f(1)(x1),  Xx(0)  =  x0, 

x2  =  ^f(0)(x2)  +  2f(1)(x2),  x2(0)  =  xi(l) 


X2(l)  = 

x(/i)  +  0(h5) 

X3(l)  : 

=  exp  (-Lf(i))  exp 

(Lf(o))  exp  (Lf(1))  (x0) 

Xl 

=  -f(1)(xi), 

xi(0)  =  x0, 

X2 

=  f(0)(x2), 

x2(0)  =  xi(l), 

X3 

=  f(1)(x3), 

x3(0)  =  x2(l) 

x3(l)  =  X 

(h)+0(h5) 
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Illustrative  Example: 


X  —  fo(t)  +  fi(t)x  +  f2(t)x2  +  f3(t)x3 


Evaluate  the  constants 


~h  Jo  (*  -  I)  mit • 


i  —  0, ...  ,3  (using  a  fourth-order  quadrature) 


Finally,  we  must  solve  the  autonomous  eqs. 


O  o 

x\  —  ao  +  +  <23X1,  ^i(O)  =  xq 

O  q 

X2  —  H-  ^1^2  H”  “1“  3^2 ?  #2(0)  —  #l(l) 


Solution  from  the  method 
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Optimization 


It  is  possible  to  improve  the  accuracy  follow¬ 
ing  different  procedures 


Suppose  the  constant  part  of  the  vector-field, 
bi,  is  the  dominant  term,  then 


1  i  bi's 

Q  =  61  —  —[6i,62]  +  *  •  *  +  cj[ 61,  •••,  61, 62]  +  •  •  • 

H - 63 - [61, 64]  -\ - [62, 63]  H - 

12  80  240 

where  the  coefficients  c*.  are  given  by 

^cfcxfc=1~i~^T=g(x). 


k>0 


x 


The  dominating  error  term  is  then 


■^1  Cj  [6l  5  •  •  •  5  61 , 62]  Cj_ j_l  [61  , 5  6l  ,  61  ,  62] 

We  can  use  additional  exponentials  to  cancel 
the  lowest  order  terms  of  E\  or  to  remove  its 
first  singularities 
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Separable  Problem: 


x  =  C(t)x  +  D(t)x 

with  [C(t),D(t)]  —  0.  Then,  we  can  consider 
the  graded  free  Lie  algebra  generated  by 


{ci,  di,  c2,  d2,  c3,  d3}  with  —  cz+di ,  i  =  1,  2,  3 
and  [ci,cy]  =  [c^,  c/j]  =  0. 

Q  —  ci  +  di  +  Y2°3  12^3_  12^Cl’  ^2]  ~  12 I^i’  C21 

[C2,  ^3]  +  fe,  C3] 


240 
1 


240 


I ’ocn([Cl,Cl^3]  T  [ci,di,C3]  +  [d1.C1.d3]  +  [dl,dl,C3]) 


360 

1 

240 
1 


([C2,  Cl,  d2]  +  [c2,  di,  C2]  +  [d2,  Cl,  d2]  +  M2,  di,  C2]) 


"^720  ([c1, Cl’  Cl’  ^2]  +  Mi,  ci,  di,  c2]  +  [ci,di,ci,d2] 

+  [ci,  di,  di,  C2]  +  [di,  ci,  ci,  d2]  +  [di,  ci,  di,  C2] 
+  [di,  di,  ci,  d2]  +  [di,  di,  di,  C2])  +  0(h7). 

exp  (si)  can  be  approximated  by 


m 


exp  (a^ici  +  QLi,2c2  +  «z,3C3)  exp  (ft,idi  +  Pi,2d2  +  ft,3ds) 


z=l 


Fourth-order:  easy 
Sixth-order:  very  complicate 
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For  the  particular  case 

x  =  Cx  +  D(t)x 
the  problem  simplifies 

£2  =  ci  +  di  +  — —  c?3  —  — -[ci,  cfe]  +  37 — ([ci,  ci,  cfe]  +  [di,  ci,  cfe]) 

12  12  360 

_ 240  (I^2’  Cl’  ^21)  ^^(tCl’Cl’Cl’d2l  +  [ci,d1,c1,d2] 
~\~[di,  ci,  ci,  cfe]  +  [rfij  dij  ci,  +  0(h7). 

and  sixth-order  methods  can  be  obtained. 

Other  cases  can  also  be  considered. 

On  the  other  hand,  for  the  autonomous  case 

x  =  Cx  +  Dx 

many  efficient  splitting  methods  are  known 

m 

eC+D  =  -Q  £aic  e(3iD 

i—  1 

Then,  for  our  composition 

m 

exp  (c^ici  +  QJz,2C2  +  Q!i,3C3)  6XP  (ft,ldl  +  ft,2^2  +  ft,3Gfe) 


we  can  take  a^i  =  a^,  f3^i  =  fa  which  solve 
many  (the  most  complicate)  order  conditions 
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Hamiltonian  Systems 


This  is  just  a  particular,  but  very  important, 
case.  Let  us  consider  i/(q,  p,£)  :  R21  xR-^M, 
where  q,  p  e  Rl,  the  Hamilton  eqs.  are 


dq  _  dH 
dt  dp  ' 
or,  equivalently 


dp 

dt 


dH 

dq 


dH 

dq 
dH 
dp 

If  we  denote  x  =  (q,  p),  it  corresponds  to  the 

dH 

particular  case  f  =  —  J — 

0x 


Lie  bracket  Poisson  bracket 


The  Hamiltonian 


H  =  T(p,t)  +  V(q,t) 


appears  in  many  problem:  is  time-dependent 
+  separable 

The  Poisson  brackets  destroy  the  separability 
Use  of  composition  Magnus  integrators 
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NUMERICAL  EXAMPLES: 
HAMILTONIAN  SYSTEMS 

Perturbed  oscillator  by  a  plane  wave 


1 

2 


0 V 2  +  q2)  +  e  (cos(q)gi(i)  +  S\n(q)g2(i)) 


with 

k 

9l(t)  =  Y  COS (wkt) 
1=1 

(  <70  =  0 

<  u>i  =  itOQ 

{  k  =  10 


k 

92(f)  =  Y  s'm(wit) 

1=1 

P0  =  11,2075 
cjq  =  1/10 
£  =  0,25  -  1,25 


We  are  interested  in  the  Poincare  section  and 
consider  different  methods.  We  start  with  a 
large  time-step,  h  and  repeat  the  computa¬ 
tions  reducing  h  until  we  get  the  correct  pic¬ 
ture.  For  this  time-step  we  measure  the  com¬ 
putational  cost 


23 


£  =  0,25 

£  =  1,25 

CPU 

N 

CPU 

N 

2EXq 

4.00 

38 

7.70 

74 

3EXq 

4.90 

58 

10.6 

121 

S * 

8.50 

38 

15.6 

71 

Srkn 

10.7 

48 

14.8 

68 

RK4 

12.0 

152 

26.5 

331 

Given  the  time  step  h  =  2tt/N  we  show  the  minimum 
value  of  N  such  that  5  <  1CT3.  For  these  values  we 
calculated  the  CPU  time  in  seconds. 
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The  Duffing  Problem 


q=  A(t)q-VqV(q,t) 
or,  equivalently 

M  =  A(t)M 

±  f  Q  \  =  (  0  V,V(q,t)  1 

dt  \  P  J  V  -M-1  0  J  \  p  / 

The  one-dimensional  equation 

q  =  eq  +  q  —  q  +£  COS  (wt) 
can  be  obtained  from  the  Hamiltonian 


5cos(cut)q 


which  is  separable  in  two  time-dependent  parts, 
being  both  of  them  solvable.  The  evaluation 
of  the  time-dependent  functions  is  the  most 
computational  costly  part 
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< 


V 


q0  —  1,75 
e  =  1 
e  10000 


P0  —  0 

e  _  1 

0  -  1000 


Errors  in  positions  for  different  splitting  meth¬ 
ods  and  the  two-exponential  CMI.  Time  step 
chosen  such  that  all  methods  require  the  same 
computational  cost. 
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The  Schrodinger  Equation 


i — u(x , 
dt  V 


/  Id2 
\  2/xdx2 


+  V(x,  t ) 


u(x ,  t) 


Spatial  semidiscretisation 


t) 


/ 


>  u(t)  =  < 


u(xo,  t) 
u(x  1,  t) 


> . 


u(xn,  t) 


Then,  we  have  to  solve 


ZU£  =  I(t)u 

with  e  CNxN  hermitic  (usually  real  and 
symmetric).  If  we  consider 

u  =  q  +  ip 

then 

q  \  =  (  0  H(t)  \  /  q  1 

dt  \  p  J  v  _HW  0  /  1  p  J 


V(x,  t)  =  D(1  -  e  ax)2  +  xf(t) 

i)  f(t)  —  Gcos(wt)  Laser  field 

ii)  fit)  —  — ,  .  Collision  with  an  atom 

J  J  v  J  cosh 2{F^wt) 
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Conclusions 


■  The  time  and  the  coordinates  in  the  vector- 
field  f(x,  f)  play  different  roles  on  the  evo¬ 
lution  of  the  system  for  many  problems 
it  is  convenient  to  treat  them  differ¬ 
ently 


■  We  have  generalized  the  Magnus  integra¬ 
tors  for  linear  systems  to  be  used  in  non¬ 
linear  problems 


■  Composition  Magnus  integrators  can  be 
used  in  tandem  with  other  geometric  in¬ 
tegrators  (of  different  order).  An  impor¬ 
tant  case  being  the  splitting  methods  for 
separable  systems.  The  final  methods  are 
still  Geometric  Integrators 


■  Additional  stages  can  be  considered  for 
optimization  purposes 
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Work  in  Progress 


■  To  analyze  the  efficiency  of  fourth-order 
methods  optimized  following  different  cri- 
terions 


■  To  build  different  sixth-order  composition 
Magnus  integrators  and  to  analyze  their 
performances 


■  To  consider  particular  cases,  like  the  sep¬ 
arable  problem,  and  to  build  efficient  com¬ 
position  methods  for  them 


■  To  look  for  interesting  problems  where 
these  methods  can  be  of  interest  (polyno¬ 
mial  vector-fields,  linear  non-homogeneous 
systems,  Riccati  equation,  separable  Hamil¬ 
tonian  systems,  some  oscillatory  problems, 
etc.) 
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